COMMENT ⊗   VALID 00018 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00003 00002	This is the stuff that does higher-level graphics, to a certain extent
C00005 00003	PROCEDURE MATX(REFERENCE REAL ARRAY A,B)
C00007 00004	PROCEDURE IDENT(REFERENCE REAL ARRAY A)
C00008 00005	PROCEDURE TRAN3(REFERENCE REAL ARRAY A REAL X,Y,Z)
C00009 00006	PROCEDURE ROT3(REFERENCE REAL ARRAY A REAL SINANG,COSANG INTEGER AXIS)
C00010 00007	PROCEDURE SCAL3(REFERENCE REAL ARRAY A REAL X,Y,Z)
C00011 00008	PROCEDURE XFORM2(REFERENCE REAL X,Y REAL ARRAY XFORM)
C00012 00009	PROCEDURE XFORM3(REFERENCE REAL X,Y,Z REAL ARRAY XFORM)
C00013 00010	INTERNAL PROCEDURE CLIP2DNORM(PROCEDURE LDRAW REAL X1,Y1,X2,Y2)
C00017 00011	PROCEDURE PDRAW(PROCEDURE LDRAW REAL X1,Y1,Z1,X2,Y2,Z2)
C00022 00012	INTERNAL PROCEDURE LINE3(PROCEDURE LDRAW REAL X1,Y1,Z1,X2,Y2,Z2)
C00025 00013	This stuff is for setting up the matrix transformations.
C00027 00014	PROCEDURE NEWVTCT
C00030 00015	PROCEDURE NEWPVP
C00032 00016	INTERNAL PROCEDURE VUPORT(REAL XPOS,YPOS,XSIZE,YSIZE)
C00033 00017	INTERNAL PROCEDURE VUWIND(REAL XAT,YAT,ZAT,XTO,YTO,ZTO,VDIST,CLPDIST)
C00034 00018	END "GRAFIX"
C00037 ENDMK
C⊗;
COMMENT This is the stuff that does higher-level graphics, to a certain extent;

ENTRY CLIP2DNORM,CLIP3DNORM,LINE3,VUPORT,VUWIND;

BEGIN "GRAFIX"

REAL ARRAY VTCT[1:4,1:4], PVP[1:4,1:4];
REAL ZNEAR,D;
REAL XE,YE,ZE,XC,YC,ZC;
REAL XW,YW,XT,YT;
PROCEDURE MATX(REFERENCE REAL ARRAY A,B);
COMMENT Multiply matrix A by B and put result in A. (Doesn't check dimensionality);
BEGIN "MATX"
REAL ARRAY TMP[ARRINFO(A,1):ARRINFO(A,2),ARRINFO(B,3):ARRINFO(B,4)];
REAL SUM;
INTEGER I,J,K;

FOR I ← ARRINFO(A,1) STEP 1 UNTIL ARRINFO(A,2) DO
    FOR J ← ARRINFO(B,3) STEP 1 UNTIL ARRINFO(B,4) DO
	BEGIN
	    SUM ← 0;
	    FOR K ← ARRINFO(B,1) STEP 1 UNTIL ARRINFO(B,2) DO
		SUM ← SUM + A[I,K]*B[K,J];
	    TMP[I,J] ← SUM
	END;
ARRTRAN(A,TMP)
END "MATX";
PROCEDURE IDENT(REFERENCE REAL ARRAY A);
COMMENT Make an identity matrix (doesn't check dimensionality);
BEGIN "IDENT"
INTEGER I,J;

FOR I ← ARRINFO(A,1) STEP 1 UNTIL ARRINFO(A,2) DO
    FOR J ← ARRINFO(A,3) STEP 1 UNTIL ARRINFO(A,4) DO
	A[I,J] ← IF I = J THEN 1.0 ELSE 0.0
END "IDENT";
PROCEDURE TRAN3(REFERENCE REAL ARRAY A; REAL X,Y,Z);
COMMENT 3-d translation;
BEGIN "TRAN3"
REAL ARRAY TRN[1:4,1:4];	COMMENT Do a real matrix translation in case the
					last column is nonzero;
IDENT(TRN);
TRN[4,1]←X; TRN[4,2]←Y; TRN[4,3]←Z;
MATX(A,TRN)
END "TRAN3";
PROCEDURE ROT3(REFERENCE REAL ARRAY A; REAL SINANG,COSANG; INTEGER AXIS);
COMMENT Do a 3-d rotation of a matrix.;
BEGIN "ROT3"
REAL ARRAY ROT[1:4,1:4];

COMMENT These formulae are on p. 335 of Newman and Sproull.;
IDENT(ROT);
CASE AXIS OF BEGIN "X Y OR Z"
["X"]	BEGIN
	ROT[2,2]←COSANG; ROT[3,3]←COSANG;
	ROT[3,2]←SINANG; ROT[2,3]←-SINANG
	END;
["Y"]	BEGIN
	ROT[1,1]←COSANG; ROT[3,3]←COSANG;
	ROT[1,3]←SINANG; ROT[3,1]←-SINANG
	END;
["Z"]	BEGIN
	ROT[1,1]←COSANG; ROT[2,2]←COSANG;
	ROT[2,1]←SINANG; ROT[1,2]←-SINANG
	END
END "X Y OR Z";
MATX(A,ROT)
END "ROT3";
PROCEDURE SCAL3(REFERENCE REAL ARRAY A; REAL X,Y,Z);
COMMENT Do a 3-d scaling of the transformation array A;
BEGIN "SCAL3"
REAL ARRAY SCL[1:4,1:4];
IDENT(SCL);
SCL[1,1]←X; SCL[2,2]←Y; SCL[3,3]←Z;
MATX(A,SCL)
END "SCAL3";
PROCEDURE XFORM2(REFERENCE REAL X,Y; REAL ARRAY XFORM);
COMMENT Do a 2-dimensional homogeneous coordinate transform;
    BEGIN "XFORM2"
	REAL XNEW,YNEW,SCALE;
	XNEW ← X*XFORM[1,1] + Y*XFORM[2,1] + XFORM[3,1];
	YNEW ← X*XFORM[1,2] + Y*XFORM[2,2] + XFORM[3,2];
	SCALE ← X*XFORM[1,3] + Y*XFORM[2,3] + XFORM[3,3];
	X ← XNEW/SCALE; Y ← YNEW/SCALE
    END "XFORM2";
PROCEDURE XFORM3(REFERENCE REAL X,Y,Z; REAL ARRAY XFORM);
COMMENT Do a 3-dimensional homogeneous coordinate transform;
    BEGIN "XFORM3"
	REAL XNEW,YNEW,ZNEW,SCALE;
	XNEW ← X*XFORM[1,1] + Y*XFORM[2,1] + Z*XFORM[3,1] + XFORM[4,1];
	YNEW ← X*XFORM[1,2] + Y*XFORM[2,2] + Z*XFORM[3,2] + XFORM[4,2];
	ZNEW ← X*XFORM[1,3] + Y*XFORM[2,3] + Z*XFORM[3,3] + XFORM[4,3];
	SCALE ← X*XFORM[1,4] + Y*XFORM[2,4] + Z*XFORM[3,4] + XFORM[4,4];
	X ← XNEW/SCALE; Y ← YNEW/SCALE; Z ← ZNEW/SCALE
    END "XFORM3";
INTERNAL PROCEDURE CLIP2DNORM(PROCEDURE LDRAW; REAL X1,Y1,X2,Y2);
COMMENT WARNING**********LDRAW procedure must have reference params.;
COMMENT This procedure clips a two-dimensional line to a square box...assuming
	that coordinates are normalized.  That is, everything is clipped to a
	box whose lower left corner is (0,0) and upper right corner is (1,1).;
COMMENT This is Cohen and Sutherland's algorithm, described in sec. 5-1 of
	Newman and Sproull's Principles of Interactive Computer Graphics, 2nd ed.;
BEGIN "CLIP2DNORM"

DEFINE  RIGHT="'2",
	LEFT="'1",
	ABOVE="'10",
	BELOW="'4",
	INSIDE="'0";

INTEGER PROCEDURE CODE(REAL X,Y);
COMMENT Returns a special code indicating the position of a point.;
BEGIN "CODE"
	INTEGER C;
	C ← INSIDE;
	IF X > 1.0 THEN C ← C LOR RIGHT ELSE IF X < 0.0 THEN C ← C LOR LEFT;
	IF Y > 1.0 THEN C ← C LOR ABOVE ELSE IF Y < 0.0 THEN C ← C LOR BELOW;
	RETURN(C)
END "CODE";

PROCEDURE HACK(REFERENCE REAL X,Y; INTEGER C);
COMMENT Clip away a protruding part of a line;
BEGIN "HACK"
COMMENT	REAL SLOPE;
COMMENT	SLOPE ← (Y2-Y1)/(X2-X1);
COMMENT Don't try to optimize common subexpressions...blows up on divide by zero.;
	IF (C LAND LEFT) THEN
	    BEGIN
		Y ← Y1 - ((Y2-Y1)/(X2-X1))*X1;
		X ← 0.0
	    END
	ELSE IF (C LAND RIGHT) THEN
	    BEGIN
		Y ← Y1 + ((Y2-Y1)/(X2-X1))*(1.0-X1);
		X ← 1.0
	    END;
	IF (C LAND BELOW) THEN
	    BEGIN
		X ← X1 - ((X2-X1)/(Y2-Y1))*Y1;
		Y ← 0.0
	    END
	ELSE IF (C LAND ABOVE) THEN
	    BEGIN
		X ← X1 + ((X2-X1)/(Y2-Y1))*(1.0-Y1);
		Y ← 1.0
	    END
END "HACK";

INTEGER C1,C2;
C1 ← CODE(X1,Y1);  C2 ← CODE(X2,Y2);
WHILE (C1 LOR C2) ≠ INSIDE DO
    BEGIN "CLIP AWAY"
	IF (C1 LAND C2) ≠ INSIDE THEN RETURN;	COMMENT Trivially invisible;
	IF C1 ≠ INSIDE THEN HACK(X1,Y1,C1);
	IF C2 ≠ INSIDE THEN HACK(X2,Y2,C2);
	C1 ← CODE(X1,Y1);  C2 ← CODE(X2,Y2)
    END "CLIP AWAY";
LDRAW(X1,Y1,X2,Y2)
END "CLIP2DNORM";
PROCEDURE PDRAW(PROCEDURE LDRAW; REAL X1,Y1,Z1,X2,Y2,Z2);
COMMENT Reference parameters because of the way SAIL works (see CLIP3DNORM);
COMMENT Make the perspective transformation, then the viewport transformation
	and the final clipping.;
BEGIN "PDRAW"
XFORM3(X1,Y1,Z1,PVP);
XFORM3(X2,Y2,Z2,PVP);
CLIP2DNORM(LDRAW,X1,Y1,X2,Y2)
END "PDRAW";

INTERNAL PROCEDURE CLIP3DNORM(PROCEDURE LDRAW; REAL X1,Y1,Z1,X2,Y2,Z2,ZNEAR);
COMMENT WARNING**********DRAW procedure must have reference params.;
COMMENT This procedure clips a three-dimensional line to a view pyramid...assuming
	that coordinates are normalized.  That is, everything is clipped to a
	pyramid bounded by the planes x=z, x=-z, y=z, y=-z, and z=ZNEAR.;
COMMENT This is Cohen and Sutherland's algorithm, described in sec. 22-6 of
	Newman and Sproull's Principles of Interactive Computer Graphics, 2nd ed.;
BEGIN "CLIP3DNORM"

DEFINE  RIGHT="'2",
	LEFT="'1",
	ABOVE="'10",
	BELOW="'4",
	INSIDE="'0",
	FRONT="'20";

INTEGER PROCEDURE CODE(REAL X,Y,Z);
COMMENT Returns a special code indicating the position of a point.;
BEGIN "CODE"
	INTEGER C;
	C ← INSIDE;
	IF X > Z THEN C ← C LOR RIGHT ELSE IF X < -Z THEN C ← C LOR LEFT;
	IF Y > Z THEN C ← C LOR ABOVE ELSE IF Y < -Z THEN C ← C LOR BELOW;
	IF Z < ZNEAR THEN C ← C LOR FRONT;
	RETURN(C)
END "CODE";

PROCEDURE HACK(REFERENCE REAL X,Y,Z; INTEGER C);
COMMENT Clip away a protruding part of a line;
BEGIN "HACK"
	REAL T;

	COMMENT This code is easier to understand if you know that the line is
		being represented parametrically as
			x = T*(X2-X1) + X1
			y = T*(Y2-Y1) + Y1
			z = T*(Z2-Z1) + Z1
		where T varies from 0 at point 1 to 1 at point 2.;
	comment the formulae were all copied from the book, one of them was verified
		using macsyma.  The z-clipping formula is not in the book and was
		derived by hand;

	IF (C LAND LEFT) THEN
	    BEGIN
		T ← (Z1+X1)/((X1-X2)-(Z2-Z1));
		Z ← T*(Z2-Z1) + Z1;
		X ← -Z;
		Y ← T*(Y2-Y1) + Y1
	    END
	ELSE IF (C LAND RIGHT) THEN
	    BEGIN
		T ← (Z1-X1)/((X2-X1)-(Z2-Z1));
		Z ← T*(Z2-Z1) + Z1;
		X ← Z;
		Y ← T*(Y2-Y1) + Y1
	    END;
	IF (C LAND BELOW) THEN
	    BEGIN
		T ← (Z1+Y1)/((Y1-Y2)-(Z2-Z1));
		Z ← T*(Z2-Z1) + Z1;
		X ← T*(X2-X1) + X1;
		Y ← -Z
	    END
	ELSE IF (C LAND ABOVE) THEN
	    BEGIN
		T ← (Z1-Y1)/((Y2-Y1)-(Z2-Z1));
		Z ← T*(Z2-Z1) + Z1;
		X ← T*(X2-X1) + X1;
		Y ← Z
	    END;
	IF (C LAND FRONT) THEN
	    BEGIN
		T ← (ZNEAR-Z1)/(Z2-Z1);
		X ← T*(X2-X1) + X1;
		Y ← T*(Y2-Y1) + Y1;
		Z ← ZNEAR
	    END
END "HACK";

INTEGER C1,C2;
C1 ← CODE(X1,Y1,Z1);  C2 ← CODE(X2,Y2,Z2);
WHILE (C1 LOR C2) ≠ INSIDE DO
    BEGIN "CLIP AWAY"
	IF (C1 LAND C2) ≠ INSIDE THEN RETURN;	COMMENT Trivially invisible;
	IF C1 ≠ INSIDE THEN HACK(X1,Y1,Z1,C1);
	IF C2 ≠ INSIDE THEN HACK(X2,Y2,Z2,C2);
	C1 ← CODE(X1,Y1,Z1);  C2 ← CODE(X2,Y2,Z2)
    END "CLIP AWAY";
PDRAW(LDRAW,X1,Y1,Z1,X2,Y2,Z2)
END "CLIP3DNORM";
INTERNAL PROCEDURE LINE3(PROCEDURE LDRAW; REAL X1,Y1,Z1,X2,Y2,Z2);
COMMENT Draw a line in 3-d.;
BEGIN "LINE3"
COMMENT Transform the points to viewing space, then to clipping space, then
	clip them against the viewing pyramid, then take the perspective
	transformation, then the viewport transformation and finally the
	screen clipping!;
XFORM3(X1,Y1,Z1,VTCT);
XFORM3(X2,Y2,Z2,VTCT);
CLIP3DNORM(LDRAW,X1,Y1,Z1,X2,Y2,Z2,ZNEAR)

END "LINE3";
COMMENT This stuff is for setting up the matrix transformations.;

PROCEDURE NEWVTCT;
COMMENT Make a new view transform/clip transform matrix.;
BEGIN "NEWVTCT"
REAL ARRAY CLP[1:4,1:4];
REAL DX,DY,DZ,DXY,DXY2,DXYZ;
INTEGER I;
IDENT(VTCT);
TRAN3(VTCT,-XE,-YE,-ZE);	COMMENT Move us to (Xe,Ye,Ze) by moving the world
					in the other direction.;
COMMENT Find description of rotations needed.;
DX ← XE-XC; DY ← YE-YC; DZ ← ZE-ZC;
DXY2 ← DX*DX+DY*DY; DXY ← SQRT(DXY2);
DXYZ ← SQRT(DXY2+DZ*DZ);

COMMENT Bring the x axis into place.;
ROT3(VTCT,DX/DXY,-DY/DXY,"Z");

COMMENT And now the y and z axes.;
ROT3(VTCT,-DZ/DXYZ,DXY/DXYZ,"X");

COMMENT Then change to the left-handed viewing coordinate system.;
FOR I ← 1 STEP 1 UNTIL 4 DO VTCT[I,2] ↔ VTCT[I,3];

COMMENT Now add on the pre-clipping transformation.;
IDENT(CLP);
CLP[1,1] ← 2*D;
CLP[2,2] ← 2*D;
MATX(VTCT,CLP)

END "NEWVTCT";
PROCEDURE NEWPVP;
COMMENT Make a new perspective and viewport transformation.;
BEGIN "NEWPVP"

COMMENT Now do the perspective transformation.;
IDENT(PVP);
PVP[3,3] ← 0; PVP[3,4] ← 2;	COMMENT 2 fixes up window size;
PVP[4,3] ← -1.0; PVP[4,4] ← 0.0;COMMENT -1 keeps z's in order;

COMMENT Move to a 0-1 coordinate space.;
TRAN3(PVP,0.5,0.5,0.0);

COMMENT Do the scaling for the viewport transformation.;
SCAL3(PVP,XW,YW,1.0);

COMMENT And finally, position the viewport.;
TRAN3(PVP,XT,YT,0.0)

END "NEWPVP";
INTERNAL PROCEDURE VUPORT(REAL XPOS,YPOS,XSIZE,YSIZE);
BEGIN
XW ← XSIZE; YW ← YSIZE; XT ← XPOS; YT ← YPOS;
NEWPVP
END;
INTERNAL PROCEDURE VUWIND(REAL XAT,YAT,ZAT,XTO,YTO,ZTO,VDIST,CLPDIST);
BEGIN
XE ← XAT; YE ← YAT; ZE ← ZAT;
XC ← XTO; YC ← YTO; ZC ← ZTO;
D ← VDIST;
ZNEAR ← CLPDIST;
NEWVTCT
END;

END "GRAFIX"